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Summary 

Finite element structural analysis based on the original 
displacement (stiffness) method has been researched and 
developed for over three decades. Although today it dominates 
the scene in terms of routine engineering use, the stiffness 
method does suffer from certain deficiencies. Various alternate 
analysis methods, commonly referred to as the mixed and 
hybrid methods, have been promoted in an attempt to com- 
pensate for some of these limitations. In recent years two new 
methods for finite element analyses of structures, within the 
framework of the original force method concept, have been 
introduced. These are termed the “integrated force method” 
and the “dual integrated force method." T . 

A comparative study was carried out to determine the 
accuracy of finite element analyses based on the stiffness 
method, a mixed method, and the new integrated force and 
dual integrated force methods. The numerical results were 
obtained with the following software: MSC/NASTRAN and 
ASKA for the stiffness method; an MHOST implementation 
for a mixed method; and GIFT for the integrated force 
methods. For the cases considered, the results indicate that, 
on an overall basis, the stiffness and mixed methods present 
some limitations. The stiffness method typically requires a 
large number of elements in the mode] to achieve acceptable 
accuracy. The MHOST mixed method tends to achieve a 
higher level of accuracy for coarse models than does the 
stiffness method, as implemented by MSC/NASTRAN and 
ASKA. The two integrated force methods, which bestow 
simultaneous emphasis on stress equilibrium and strain com- 
patibility, yield accurate solutions with fewer elements. in a 
model. The full potential of these new integrated force methods 
remains largely unexploited, and they hold the promise of 
spawning new finite element structural analysis tools. 

Introduction and Overview 

The field of finite element analysis for structures, based on 
the original stiffness method and the more contemporary mixed 
and hybrid methods, has made great strides during the past 
three decades. General purpose finite element software such 
as MSC/NASTRAN (ref. 1) and ASKA (ref. 2), based on the 
stiffness method, and MHOST (ref. 3), based on a mixed- 
iterative method, are examples of structural analysis tools 
available today. The current generation of finite element 


analysis software, coupled with modern computer hardware, 
provides the capability to solve challenging engineering 
problems that require extensive numerical calculations. Despite 
their popularity and prominence, the current finite element 
analysis methods are not free from deficiencies, and oppor- 
tunities for improvement appear to exist. In an attempt to 
compensate for some of the limitations, two new formulations 
within the framework of force method concepts have been 
introduced during the past few years. These are termed the 
“integrated force method"" (IFM) and the “dual integrated 
force method"’ (IFMD). This report examines the accuracy 
of finite element structural analysis via the IFM and IFMD 
versus analysis by the stiffness and the mixed and hybrid 
methods. 

An overall qualitative assessment of the various analysis 
methods can be attempted from a consideration of the universal 
equilibrium equations, which represent the force or stress 
balance conditions. The force equilibrium conditions, in the 
general context of finite element analysis, give rise to an 
unsolvable indeterminate system of equations with a greater 
number of unknown forces than the number of such equations. 
The equilibrium equations, being indeterminate, cannot be 
solved for the unknown forces, except for the trivial statically 
determinate case. Because of the indeterminancy, various 
alternative methods have been devised for stress analysis of 
indeterminate structures. The methods available for finite 
element analysis that are of interest in this study are briefly 
summarized in the next section. The nomenclature for the 
analysis method adapted in this paper is based on the primary 
unknown of the formulation; these unknowns are defined in 
table l and illustrated in figures 1 and 2. 

The Integrated Force Method— A Direct Force Method 

In the direct force method all of the internal forces are treated 
as the primary unknowns and are directly computed by solving 
a_ set of simultaneous equations, A solvable system of equations 
is obtained by augmenting the rectangular system of 
equilibrium equations with another rectangular system of 
equations expressed in terms of the same unknown forces. The 
augmenting system represents the strain compatibility condi- 
tions. The total system resulting from the concatenation of the 
force equilibrium equations and strain compatability conditions 
is a solvable set of n equations in n unknowns, the solution 
of which directly yields all n internal forces. This direct force 
method, which bestows simultaneous emphasis on both the 
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TABLE I. -METHODS OF STRUCTURAL ANALYSTS AND ASSOCIATED 
VARIATIONAL FUNCTIONALS 


Name of method 

Primary variables 

Variational 

functional 

Elasticity 

Structures 

Elasticity 

Structures 

Completed Beltrami- 
Michell formulation 

Integrated force 
method (IFM) 

Stresses 

Forces 

IFM variational 
functional 

Navier formulation 

Stiffness method 

Displacements 

Deflections 

Potential energy 

Airy formulation 

Classical force 
method 

Stress function 

Redundants 

Complementary 

energy 

Mixed formulation 

Reissner method 

Stresses and 
displacements 

Forces and 
deflections 

Reissner 

functional 

Total formulation 

Washizu method 

Stresses, strains, 
and displacements 

Forces, 
deformations, 
and deflections 

Washizu 

functional 



Figure L— Force and displacement methods. 


equilibrium equations and the compatibility conditions and 
solves for the forces directly, is the integrated force method 
(TFM) (refs. 4 to 19). The additional key ingredient for the 
IFM, which parallels the completed Bel trami-Michell formula- 
tion of elasticity (refs. 5, 14, 17, and 20), is the explicit 
formulation of the global strain compatibility conditions of 
finite element models. These compatability conditions of finite 
element models, which are analogous to St. Venant’s strain 
formulation of elasticity, have been divided into interface, 
cluster, and boundary compatibility conditions (refs. 8, 10, 
and 11). They enforce deformation balance (1) along the 
element interface ’ (2) for a cluster of elements, and (3) along 



Figure 2.— Hybrid and mixed methods. 


the constrained segment of the boundary of the discrete model. 
The IFM was not developed in the formative 1960’s because 
the concept of the compatibility conditions augmenting the 
equilibrium equations for indeterminate structures had not yet 
been recognized. 

Redundant Force Method 

Despite the nonavailability of an explicit, computer- 
automated compatibility formulation, a second classical analysis 
method, known as the redundant force method (refs. 21 and 
22), was developed. The redundant force method was first 
formulated by Maxwell in the mid 1800's and remained the 
analysis method of choice for about a century. Recognizing 
the indeterminate nature of the equilibrium equations, Maxwell 
introduced the ingenious concept of redundant forces and their 
use in analysis. In redundant analysis, elements are “cut” to 
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create a determinate system; then ad hoc compatibility is 
restored by closing the “gaps”. This key procedure yields a 
system of r equations in terms of r unknown redundant forces. 
Redundant analysis yields redundant forces that are treated as 
external loads on the auxiliary determinate structure. The 
indeterminate analysis is completed by invoking the principle 
of superposition on the determinant structure. 

The redundant force method requires the selection of an 
auxiliary determinate basis structure and corresponding 
redundant forces, which is the major difficulty of this classical 
analysis method. Prior to the availability of computers and 
the computerization of matrix methods, redundants were 
identified manually; such a process depended on the subjective 
experience and judgment of the structural analyst. Sub- 
sequently, redundant identification was automated, at least in 
principle (refs. 23 to 31), on the basis of linear algebra 
concepts such as rank, column combination and diagonal 
dominance of the equilibrium matrix, and a self-equilibrating 
stress state. Such concepts, although analytically elegant, 
lacked the physical features of the compatibility conditions of 
finite element models (i.e., deformation balance among 
element interfaces, clusters, and constrained boundary 
segments) and the desirable numerical features, such as 
bandwidth and conditioning of the compatibility matrix. As 
a result, the redundant force method failed for structures of 
any complexity. Despite the attention of and earnest efforts 
by prominent researchers (refs. 23, 25, 27, and 30), the Air 
Force, and NASA, the redundant force method never became 
an integral part of the well-known finite element software 
NASTRAN. The original intent of NASA was to provide for 
both force and stiffness methods in NASTRAN; however, only 
the stiffness method implementation now exists. 

To illustrate the complexity associated with automatic redun- 
dant selection, consider as an example a plate flexure problem 
(see appendix) in which the plate is discretized by two finite 
elements and the model has m = 12 displacement unknowns 
and n = 18 force unknowns. The maximum possible number 
of redundant force systems is given by 

n\ 


or c r nVdX =49,504 in this case, of which probably only one 
is the desired canonical set. The maximum possible number 
of redundant systems from which a canonical set can be 
selected increases rapidly for more complex structures. For 
example, for m- 15 displacement unknowns and n = 25 force 
unknowns, the maximum number of possible sets exceeds 
3 million. Attempts have been made to reduce such large 
numbers of choices for redundants; however, the problem has 
not been satisfactorily resolved because of its intrinsically 
difficult nature. Overall, the classical redundant force method 
as a computerized method of analysis outlived its usefulness 
and was abandoned during the early stages of the development 


of computerized structural analysis technology. Its inclusion 
in the discussion here is for completeness, but it is not 
considered further in this study. 

Stiffness Method 

The statically indeterminate nature of the equilibrium 
equations and the nonexistence of a strain compatibility formu- 
lation led to the direct displacement analysis, or original 
stiffness, method (ref. 21). In the stiffness method, first 
formulated by Clebsch (1833 to 1872), the equilibrium equa- 
tions are written in terms of displacements that when 
augmented with the displacement continuity conditions give 
rise to an adequate system of equations with which to solve 
for the unknown displacements (see the appendix). From the 
known displacement state, the forces and stress parameters 
are obtained as secondary variables by differentiation or its 
equivalent, which could tend to degrade the accuracy of the 
stress predictions. The stiffness method, which generally 
requires extensive computations, was not popular until the 
emergence of high-speed computers. Since the dawn of the 
computer age (for the past three decades), the stiffness method 
has been extensively researched throughout the world, and 
today it dominates the engineering analysis scene. 

Mixed and Hybrid Methods 

Because of the limitations of the stiffness method, especially 
with respect to the accuracy of stress solutions, two other 
approaches have been devised for finite element analysis of 
structures. The method that considers both stresses and 
displacements as simultaneous unknowns is referred to as the 
hybrid method (refs. 31 to 33 and fig. 2). The other method, 
known as the mixed method (refs. 3 and 34 to 38 and fig. 2), 
treats displacements, stresses, and strains as simultaneous 
unknowns. Neither the hybrid nor the mixed method imposes 
strain compatibility conditions explicitly. Rather, these two 
methods systematically combine the equilibrium equations, the 
displacement continuity conditions, the kinematic relations, 
and the material constitutive relations in one form or another. 
The hybrid and the mixed methods, although generally more 
computationally demanding than the direct methods (for 
comparable discretizations), do not, however, contain any 
fundamentally new ingredient that is not already present in 
the force and displacement methods. 

The various analysis methods that have been associated with 
underlying variational functionals are summarized in table I 
and depicted in figures 1 and 2. 

In these figures the following relationships are shown: 

(1) The displacement method explicitly utilizes the equi- 
librium equations written in terms of displacements, which 
are augmented with the displacement continuity conditions 
(fig. 1). 

(2) The hybrid method includes equilibrium equations, stress 
displacement relations, and the displacement continuity 
conditions (fig. 2). 
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(3) The mixed method uses equilibrium equations, strain 
displacement relations, the constitutive law, and the displace- 
ment continuity conditions. (Note that the equilibrium 
equations are considered in all the methods.) 

(4) The IFM is the only formulation that makes use of the 
strain compatibility conditions along the interface, field, and 
boundaries of the finite element model in addition to the 
equilibrium equations (fig. 1). 

(5) The displacement method and the hybrid and mixed 
methods do not explicitly make use of the strain compatibility 
conditions (figs. I and 2). 

If variables are classified with respect to the universal law 
of equilibrium, then forces are its primal variables, and 
displacements are its dual variables. On this basis, the IFM 
is the primal analysis method since its unknowns are the forces, 
and the stiffness method is the dual method since its unknowns 
are the displacements. At present, the dual displacement, or 
stiffness, method has been exhaustively researched, and its 
potential has been exploited to the extent that the method may 
have reached a plateau in its development. Conversely, the 
primal analysis method, with the emergence of the two new 
integrated force methods, appears to hold considerable 
potential for further development. 

Finite Element Analysis 

In the discrete finite element analysis technique, the element 
characteristics and external loads are lumped at the nodal points 
of the model, and the governing equations are written with 
respect to these grid points. The solution obtained by finite 
element analysis should satisfy the two fundamental axioms of 
structural mechanics (i.e., the satisfaction of the force equilib- 
rium equations and the compliance of the strain compatibility 
conditions), at least with reference to the grid points of the 
model. Even though the stiffness method depends heavily on 
the state of equilibrium at the nodal points, it is commonly 
observed that at those very cardinal points stresses recovered 
from the nodal displacements often violate equilibrium. The 
mixed method of MHOST attempts to compensate for this 
limitation through an iterative solution process that equilibrates 
stresses at the node points. Although stress equilibrium 
imbalance has been researched, the problems have not been 
resolved to complete satisfaction. The IFM, in which the internal 
force parameters are explicitly constrained to simultaneously 
satisfy both the equilibrium equations and the compatibility con- 
ditions with reference to the grid points, is an attempt to obtain 
accurate stresses even at the nodes of the finite element model . 

The purpose of this report is to examine the accuracy of the 
different methods of finite element analysis. To accomplish this, 
the results obtained for several test problems by different methods 
were scrutinized with regard to the relative performance. 

Numerical solutions for the test cases were obtained with 
the following finite element software: 

GIFT— Based on the theories of the IFM and the dual IFMD, 
GIFT is a modest program developed for research purposes. 


MSC/ NASTRAN .— This program is one of the most widely 
used stiffness-method-based codes available today. 

MHOST .— Designed especially for nonlinear analysis, this 
program provides a versatile analysis capability based on a 
mixed-iterative formulation. 

ASKA .— The ASKA program, developed in Europe, is also 
based on the stiffness method. It is used here mainly for 
numerical verification with the MSC/NASTRAN code. 

Hybrid Method .— Although obtained independent of this 
study, solutions from a hybrid method using element HMPL5 
(ref. 34) are included here for the sake for completeness. 

This report does not attempt to elaborate on the theoretical 
details of the different analysis methods. References 4 to 19 
can be examined for representative research results on the 
theory of the IFM for elastic continua, finite element analysis, 
and design optimization. 

The subject matter of this report is presented in Four sections: 
the basic equations of the analysis; the test cases and results; 
a discussion of the results; and conclusions. In addition, to 
illustrate the calculation sequence for the force and the stiffness 
methods, an example is provided in appendix A. Symbols are 
defined in appendix B. 

Basic Equations of the Methods 

This section summarizes the governing equations of the analysis 
methods investigated here, namely, (1) integrated force method, 
(2) dual integrated force method, (3) stiffness method, and 
(4) mixed- iterative method. The equations of the hybrid method 
are not presented. For detailed examination of the theories of 
the methods, references 1 to 8, 20, and 21 are suggested. 

Integrated Force Method and Dual Integrated 
Force Method 

In the integrated force method (IFM), the internal forces 
are taken as the primary unknowns and the displacements are 
obtained by a back calculation operation. The dual integrated 
force method (IFMD) is derived from the equations of the IFM 
by eliminating internal forces in favor of displacements. The 
primal variables of the IFMD are, then, the displacements from 
which forces are recovered in secondary operations. 

In the IFM a discretized structure for the purpose of analysis 
is designated by attributes («,/n), which denote the number 
of force and displacement degrees of freedoms, respectively. 
In the IFM analysis a governing set of n equations is expressed 
in terms of n unknown internal forces [F]. The system of n 
equations is obtained by augmenting the set of m force 
equilibrium equations 

[Bi :fj = [p] 

with the set of r = n — m strain compatibility conditions 
[C][G](F| = [SRj 
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as follows: 


[B] 

[C][G] 


m 



or [S][Fi = (P*j 


(1) 


where [B] is the mXn equilibrium matrix, [C] is the rxn 
compatibility matrix, [G] is the nxn concatenated flexibility 
matrix that links deformations [0] to forces [F] as 

[0) = [G][Fj 

(Fj is the 1 internal force vector, (Pj is the mX 1 external 
load vector, and [5R] is the rx 1 effective initial deformation 
vector defined by 


ators of elasticity and the matrices of the TFM. Mathematically 
speaking, the derived operators and matrices of other 
formulations can possess characteristics (i.e., numerical 
norms, spectral radii, and stability of equation systems) no 
more superior than the basic unsymmetrical operators of 
elasticity theory or matrices of the IFM, even when the derived 
operators and matrices become symmetric (refs. 6 and 7). 

The governing equations of the IFMD are generated from 
the IFM equations (1) and (2) by mapping forces into 
displacements and vice versa. The key equation of the IFMD, 
wherein nodal displacement unknowns [X] become the primary 
variables and are linked to the external loads [P], resembles 
the familiar stiffness equation and is given as 

[K,](X) = (Pj (3) 


[«RJ = " [C](j8 0 ] 

where [j3 0 ] is the nx 1 initial deformation vector, and [S] is 
the nxn IFM governing matrix. The matrices [B], [C], [G], 
and [S] are banded and have full row ranks of m, r, n , and 
n , respectively. The solution of equation (1) yields n internal 
forces. The displacements are obtained from the forces in a 
back calculation operation expressed as 

[XJ = [J]([G](F] + [&>}] (2) 

where [J] is the deformation coefficient matrix defined as the 
first mxn partition of [[S] “ '] . Equations (1) and (2) 
represent the two key relations of the IFM for finite element 
analysis that are needed to calculate the forces and 
displacements, respectively. 

In terms of fundamental operators, an analogy can be made 
between the IFM and the theory of elasticity (ref 40). The 
three fundamental operators of elasticity are (1) the equilibrium 
operator of Cauchy, which relates stresses to external loads; 
(2) the compatibility operator of St. Venant, which controls 
components of strain; and (3) the material constitutive tensor 
of Hooke, which relates strains to stresses. Likewise, the IFM 
has three operators that are equivalent to the operators of the 
elasticity theory. The operators, which become matrices in 
the context of finite element analysis, are (1) the equilibrium 
matrix [B], which links internal forces to external loads; (2) 
the compatibility matrix [C], which governs the deformations; 
and (3) the flexibility matrix [G], which relates deformations 
and forces. Both the equilibrium and the compatibility 
operators of elasticity and the corresponding matrices of the 
IFM are nonsymmetrical, whereas the material constitutive 
tensor and the flexibility matrix are symmetrical. Governing 
operators of other formulations (e.g., Navier’s displacement 
formulation, Airy’s stress function formulation, Reissner’s 
hybrid formulation, or the Hu-Washizu’s mixed formulation) 
and the matrices of other discrete analysis methods (such as 
the stiffness, redundant force, mixed, and hybrid methods) 
are, in principle, derivable from the basic unsymmetrical oper- 


where [K v ] is a matrix defined by the first mxm partition of 
the matrix product [[S][G] _ '[S] r ] and is referred to as the 
pseudo-stiffness matrix. 

For the element types that have been formulated to date 
(including rectangular membrane and flexure elements, 
triangular membrane and flexure elements, and solid brick and 
tetrahedral elements), we have observed that, for consistent 
force and displacement field assumptions, the attributes of the 
pseudo-stiffness matrix (such as symmetry, dimension, and 
sparsity) are identical to those of the conventional stiffness 
matrix [K]. Only the magnitudes of nonzero coefficients of 
the two matrices [K] and [KJ differ. 

Once displacements are obtained as the solution to equa- 
tion (3), forces can be obtained by back calculation as 

(Fj = [G 7 ](X] (4) 

where the nXm force coefficient matrix [Gy] is 
nonsymmetrical and is defined in terms of the product of the 
inverse of the concatenated flexibility matrix [G] _l and the 
first n X m partition of the transpose of the IFM governing 
matrix [S] r (defined in eq. (1)). Since the flexibility matrix 
[G] is the block diagonal concatenation of the corresponding 
element matrices, its inverse is inexpensive to compute, and 
calculating forces from displacements by using equation (4) 
requires only a small fraction of the total computations 
necessary for the entire analysis. 

Since equations (1) and (2) of the IFM are mathematically 
equivalent to equations (3) and (4) of the IFMD, the forces, 
displacements, and deformations obtained by either of the 
methods are identical; thus, 


[F]ifm “ (FJifmd 

(5a) 

PCjiFM — (X)IFMD 

(5b) 

($IFM = (0) IFMD 

(5c) 
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The relations given by equation (5) have also been observed 
numerically; that is, the numerical results obtained for each 
test case satisfied equation (5) as expected. 

The Stiffness Method 

The governing equation of the stiffness method, wherein the 
primary variables [X] (the nodal displacements of the finite 
element model) are linked to external loads [P] through the 
stiffness matrix [K], can be symbolized as 

[K]fX] = [PJ (6) 

where [K] is the symmetrical stiffness matrix of dimension 
m x m . 

Unlike the integrated force technique (eq. (2) or eq. (4)), 
the stiffness method has no single expression that can be used 
in calculating stress parameters from displacements by back 
calculations. The equivalent of differentiation and a series of 
numerical operations are required to generate the deformation 
and force variables from the nodal displacements. 


The matrix [K] is the standard stiffness matrix as in equation 
(6); [E,J is an integrated form of the discrete gradient matrix; 
[GJ is the material elastic constitutive matrix; and [C m ] is 
the strain projection matrix. The vectors (uj, [e], [a], and (Pj 
represent the nodal displacements, strains, stresses, and loads, 
respectively. The matrices [NJ and [S a ] comprise the 
interpolating polynomials (shape functions) for the strain and 
stress fields, respectively. 

The iterative process of the MHOST mixed-iterative strategy 
is as follows: 

Step 1: Initial stiffness solution 

[XJ = [K]-'(P] (9a) 

Step 2: Nodal displacement update 

[X] n+ 1 = [u] fl + [K] - 1 [[Pj - [E] r (<] (9b) 

Step 3: Nodal strain projection 


The MHOST Mixed-Iterative Method 

The MHOST finite element code (refs. 3, 36, and 37) 
implements a mixed-iterative method derived from an 
augmented Hu-Washizu variational principle, and it employs 
an equal-order interpolation of the fields, displacement, strain, 
and stress, which are represented consistently as nodal variables. 
The mixed equations of the general Hu-Washizu formulation 
are augmented with the conventional stiffness equation and 
solved indirectly; that is, the stiffness equation serves as a 
preconditioner for the iterative recovery of the mixed solution. 
This avoids the computational penalty of a direct solution of 
the mixed equation system in which all three fields are treated 
as simultaneous unknowns. The governing equations of the 
MHOST mixed-iterative method are expressed as 


{*) n+1 = [C] -l [E][Xj" +l (9c) 

Step 4: Nodal stress recovery 

M n+1 = [C]- r [G](cJ n+l (9d) 

Step 5: Evaluation of the nodal equilibrium residual 

(r] n+l = {Pj - [E] r [<r) n+, {r] n+1 = [P] - [E] r K + l (9e) 

The process iterates on steps 2 to 5 to reduce the residual vector 
defined by equation (9e) to an acceptable level. 

The initial solution given by equation (9a) is the standard 
stiffness solution, and in MHOST terminology, it is referred 
to as the MHOST/uniterated solution. The converged solution 
from steps 2 to 5 is referred to as the MHOST/iterated solution. 


‘ [K] [0] [E ,„] 7 ' 


M 


[P - Kuj 

[0] [G,„] [-CJ r 


W 

= 

(0j 

. [EJ [-CJ [01 . 


f {„)} 

1 1 

10] 


where 


[EJ = 

yN ff i r [B]jn 

(8a) 

[GJ = 

j n [N e ] r [D][N ( ]</n 

(8b) 

[CJ = 

LfN.rtN f ]dQ 

(8c) 


Numerical Results 

Numerical results for test problems obtained by different 
methods are presented in this section. The attributes of the 
finite elements used by the different software are presented 
next. 

MSC/NASTRAN Elements QUAD-4 and TRIA-3 

Four-node QUAD-4 and three-node TRIA-3 elements of 
MSC/NASTRAN were used in this study, with the displace- 
ment degrees of freedom constrained in such a way as to 
separately obtain the membrane and bending responses. For 
bending response, the degrees of freedom are restricted to a 
transverse translation and the two rotations; the QUAD-4, 
then, is a 12-degree-of-freedom element, and TRIA-3, a 
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9-degree-of-freedom element. For membrane response, the 
QUAD-4 element has eight degrees of freedom, that is, two 
in-plane translations for each of its nodes. 

ASKA Elements QUAD-4, TRIB-3, and TUBA-3 

The ASKA finite element software also has a QUAD-4 
element that was used for this study. The attributes of the 
QUAD-4 element of ASKA are identical with respect to nodes 
and degrees of freedom to those of the QUAD-4 element of 
MSC/NASTRAN, Two triangular elements of the ASKA soft- 
ware were used to examine the difference in performance of 
higher order elements in finite element calculations. The 
element TRIB-3 for flexural response has three degrees of 
freedom per node, consisting of a transverse translation and 
two rotations. Element TUBA-3 is a higher order triangular 
element, which for bending response alone has six degrees 
of freedom per node, consisting of one transverse translation, 
two rotations, and three curvatures. As will be seen, the 
increase from three to six degrees of freedom per node did 
not significantly improve the accuracy in the cases considered. 

GIFT Elements PLB4SP, MEMRSP, and PLB3SP 

The elements of GIFT software used for this study were the 
four-node plate-bending element PLB4SP, the three-node 
plate-bending element PLB3SP, and the four-node membrane 
element MEMRSP (the same element name is used for both 
the IFM and the IFMD). 

The IFM element PLB4SP has three force degrees of 
freedom per node, consisting of one shear force and two 



Nodal displacements 
Membrane element MEMRSP 

Figure 3.— Four-node 


moments, whereas the IFMD element PLB4SP has three dis- 
placement degrees of freedom per node, consisting of a 
transverse translation and two rotations. The PLB4SP element 
for IFMD corresponds to the restrained QUAD-4 elements 
of MSC/NASTRAN and ASKA. The force and displacement 
degrees of freedom of the PLB4SP elements are depicted in 
figure 3. Likewise, the restrained three-noded triangular 
element TRIA-3 of MSC/NASTRAN and the IFM/IFMD 
element PLB3SP are equivalent; TRIA-3 corresponds to 
translation along the transverse direction and rotations along 
the two in-plane axes, whereas the IFM element PLB3SP 
represents nodal forces along those directions for the IFM 
element. The MEMRSP element is a four-node rectangular 
membrane element. For the IFM the MEMRSP element con- 
tains two force degrees of freedom per node, representing the 
two membrane forces along the coordinate axes, and for the 
IFMD it has two displacement degrees of freedom per node. 

MHOST Elements SH75 and PS151 

The MHOST element SH75 used in this study is a four-node, 
bilinear, isoparametric, quadrilateral element based on 

Reissner-Mindlin plate and shell theory (refs. 3 and 38). It 
has six displacement degrees of freedom per node (three 
translations and three rotations). The element is formulated 
in terms of nine generalized deformations, consisting of strains 
and curvatures (e A , e yr e z , y xy , y vz , y xz , k x , k v , k X} ), and nine 
generalized stress resultants (N x , N yt N xy , S YZ , S KZ , M x> M y , 
M xy , * z ). 

The MHOST plane stress element PS151 is a four-node, 
bilinear, isoparametric, quadrilateral element based on 


s 4 s 7 



Nodal displacements 
Flexure element PLB4SP 

and flexure elements. 
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independent strain interpolation. The nodal variables for the 
element include two displacements (u x> m v ), three strains (e u , 
e yv , 7. n ), and three stresses (cr xy , o xy r vv ). 

Overall, the elements QUAD-4 and TRIA-3 of MSC/ 
NASTRAN, QUAD-4 and TRIB-3 of ASKA, and PLB4SP, 
MEMRSP, and PLB3SP of GIFT are ‘ ‘ordinary’ ’ elements 
with three degrees of freedom for bending response and two 
degrees of freedom for membrane response, and they can be 
considered equivalent to one another. The elements TUBA-3 
of ASKA, and SH75 and PS151 of MHOST can be considered 
high-precision elements, either because their nodal degrees 
of freedom exceed those of the normal elements or because 
an iterative residue-controlling scheme is adopted as in the 
MHOST/iterative scheme. 

The test cases considered for this study are summarized as 
follows. 

Case I— Analysis of the Cantilever Beam 

The cantilever beam, shown in figure 4, represents a typical 
finite element test problem. The beam is made of an isotropic 
material, and its parameters are as follows: 


Length, a , in 24 

Depth, d , in 2 

Thickness, /, in 0.25 

Young’s modulus, E, ksi 30 000 

Poisson’s ratio, v 0.3 

Magnitude of transverse concentrated 
load at each of two free end nodes, lb 100 


The theoretical solutions for the cantilever beam are as follows 
(ref. 42): displacement at the tip of the beam is 

5y = 0. 18432 in. (10a) 

shear force at any location along span x of the beam is 

V y = 200 lb (10b) 

and bending moment along span x of the beam is 

M x = 200(24 - jc)lb-in. (10c) 

The beam was discretized as shown in figure 4(b). It was 
analyzed by using the quadrilateral elements QUAD-4 of 
MSC/NASTRAN, element MEMRSP of GIFT, and elements 
SFI75 and PS151 of MHOST. The computed results for dis- 
placement and stress were normalized with respect to the 
theoretical solutions. The displacement and stress results along 
with the equilibrium imbalance at the nodes of the finite 
element model are presented in tables II to IX. However, nodal 
stresses obtained by stiffness methods (MSC/NASTRAN and 
ASKA) were ambiguous (ref. 39); therefore these are not 
included in table IV. The convergence of the tip displacement 
solution with respect to the number of finite elements in the 
discrete beam model is depicted in figure 5. The stiffness 
(eq. 6) and the pseudo-stiffness (eq. 3) coefficients for a 
12-element model are given in table X. 


-a = 24 in. 



P/2 


P/2 



(a) Geometry and boundary conditions. 

(b) Finite element model. 

Figure 4.— Cantilever beam analysis— Case I. 
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TABLE II. -COMPARISON OF IFM, MSC/NASTRAN, AND MHOST 
CANTILEVER BEAM NORMALIZED TIP DISPLACEMENTS 3 


Number of 
elements, 
n 

Normalized displacement 

GIFT/IFM 

MEMRSP 

GIFT/IFMD 

MEMRSP 

MSC/NASTRAN 

QUAD-4 

MHOST 

SH75 

Uniterated 

solution 

Iterated 

solution 

[ 

0.755 

0.755 

0.614 

0.678 

0.678 

2 

.942 

.942 

.858 

.855 

1.024 

3 

.977 

.977 

.889 

.888 

.955 

4 

.989 

.989 

.900 

.900 

.945 

6 

.998 

.998 

.908 

.908 

.927 

8 

1.000 

1.000 

.911 

.911 

.921 

10 

1.002 

1.002 

.912 

.912 

.919 

12 

1.002 

1.002 

.913 

.913 

.918 




.914 

.914 

.914 

10 

'SA 



.914 

.914 

.914 

j 48 



.914 

.914 

.914 


a Unil> represents analytical solution. 


TABLE m. -COMPARISON OF IFM AND MHOST 
ELEMENT PS 151 CANTILEVER BEAM 
NORMALIZED TIP DISPLACEMENTS 3 


Number of 

Normalized displacement 


elements, 

n 

GIFT/IFM 

MEMRSP 

GIFT/IFMD 

MEMRSP 

MHOST 

PS151 




Uniterated 

solution 

Iterated 

solution b 

1 

0.755 

0.755 

0.755 

0.755 

2 

.942 

.942 

.942 

.986 

3 

.977 

.977 

.977 

.989 

4 

.989 

.989 

.989 

.994 

6 

.998 

.998 

.998 

.999 

8 

1.000 

1.000 

1.001 

1.001 

10 

1.002 

1.002 

1.002 

1.002 

12 

1.002 

1.002 

1.003 

1.003 

16 

24 



1.004 

1.004 

1.004 

1.004 


a Unity represents analytical solution. 
b For one iteration. 


TABLE IV — COMPARISON OF IFM AND MHOST ELEMENT 
SH75 CANTILEVER BEAM NORMALIZED 
STRESSES AT SUPPORT 3 


Number of 
elements, 
n 

Normalized stress 

GIFT/IFM 

MEMRSP 

GIFT/IFMD 

MEMRSP 

MHOST 

SH75 

Uniterated 

solution 

Iterated 

solution 

1 

1.000 

1.000 

0.494 

0.494 

2 





.748 

.995 

3 





.832 

.943 

4 





.874 

.997 

6 





.916 

1.003 

8 





.937 

.985 

10 





.950 

.988 

12 





.958 

.990 

16 





.969 

.982 

24 


r 



.979 

.988 


a Unity represents analytical solution. 
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TABLE V. -COMPARISON OF IFM AND MHOST ELEMENT 
PS 15 1 CANTILEVER BEAM NORMALIZED 
STRESSES AT SUPPORT 3 


Number of 

Normalized stress 

elements. 






: 

n 

GIFT/IFM 

GIFT/IFM D 

MHOST 


MEMRSP 

MEMRSP 

PSI5I 






Uniterated 

Iterated 






solution 

so!ution h 

1 

1.000 

1.000 

0.500 

0.500 

2 





.750 

.838 

3 





.833 

.881 

4 





.875 

.912 

6 





.917 

.944 

8 





.938 

.960 

10 





.950 

.969 

12 





.958 

.975 

16 





.969 

.982 

24 





.979 

.988 


a Unii> represent analytical solution, 
ror one iteration. 


TABLE VI. -COMPARISON OF IFM AND MHOST ELEMENT 
SH75 CANTILEVER BEAM EQUILIBRIUM 
IMBALANCES AT POINT A a 


Number of 

Equilibrium imbalance, percent 

elements. 






— 

n 

GIFT/IFM 

GIFT/IFM D 

MHOST 


MEMRSP 

MEMRSP 

SH75 






Uniterated 

Iterated 



.. 



solution 

solution b 

2 

0 

0 

49.839 


3 





12.478 

— 

4 





8.322 

0.104 

6 





4.994 

.515 

8 





3.567 

.380 

10 





2.774 

.312 

12 





2.270 

.263 

16 





1.664 

.825 

24 





1.085 

.518 


a Unity represents analytical solution. 
b Fnr one iteration 


TABLE VIL— COMPARISON OF IFM AND MHOST ELEMENT 
PS 15 1 CANTILEVER BEAM EQUILIBRIUM 
IMBALANCES AT POINT A a 


Number of 
elements, 
n 

Equilibrium imbalance, percent 

GIFT/IFM 

MEMRSP 

GIFT/IFMD 

MEMRSP 

MHOST 

PS151 

Uniterated 

solution 

Iterated 

solution h 

2 

0 

0 

0 

-5.847 

3 







1.807 

4 







1.130 

6 







0.921 

8 







.833 

10 







.752 

12 







.682 

16 







.570 

24 







.424 










a Unity represents analytical solution. 
Vor one iteration. 


TABLE vrn. -IMPROVEMENT IN MHOST ELEMENT 
SH75 CANTILEVER BEAM SOLUTION 


Number of 
elements, 
n 

Percentage 

improvement 

Computational 

penalty 

Support 

stress 

Tip 

displacement 

Number of 
iterations 

Normalized 
extra time 

1 



0 


2 

24.70 

16.90 

2 

1.877 

3 

11.10 

6.70 

2 

1.022 

4 

12.30 

4.50 

3 

2.400 

6 

8.70 

1.90 

5 

3.491 

8 

4.80 

1.00 

2 

2.00 

10 

3.80 

.70 

2 

2.04 

12 

3.20 

.50 

2 

2.037 

16 

1.30 


1 

1.650 

24 



1 

1.000 

Theory 

100.0 

100.0 

— 

— 


(0 





TABLE IX. -IMPROVEMENT IN MHOST ELEMENT 
PS 15 1 CANTILEVER BEAM SOLUTION 


Number of 
elements, 
n 

Percentage 

improvement 

Computational 

penalty 

Support 

stress 

Tip 

displacement 

Number of 
iterations 

Normalized 
extra time 

1 


___ 

0 



2 

8.8 

4.4 

I 

1.20 

3 

4.8 

1.2 

1 

1.33 

4 

3.7 

.5 

1 

1.23 

6 

2.7 

.1 

I 

1.33 

8 

2.2 

— 

1 

1.33 

10 

1.9 

— ! 

1 

1.27 

12 

1.7 

— 

1 

1.29 

16 

1.3 

— 

1 

1.26 

24 

.9 

— 

1 

1.37 


TABLE X. -CANTILEVER BEAM STIFFNESS MATRIX 
COEFFICIENTS FOR 12-ELEMENT MODEL OF 
STIFFNESS MATRIX DIMENSION (48,48) 


Stiffness 

GIFT/IFMD 

MSC/NASTRAN 

Difference, 

matrix 

(MEMRSP) 

(QUAD-4) 

percent 

coefficients, 




K v 




(48th row) 




^48.01 t0 ^48.20 

0 

0 

0 

^48.21 

-1.339x 10 6 

-1.339X10 6 

0 

^48,22 

-2.157 

-2.095 

2.874 

^48,23 

-.103 

-.103 

0 

^48.24 

-1.964 

-2.026 

-3.157 

K48. 25 t0 ^48 .44 

0 

0 

0 

^48,45 

.134 

.103 

23.134 

* 

OC 

& 

.714 

.652 

8.863 

^48,47 

1.339 

-1.339 

0 

^48 .48 

3.407 

3.468 

-1.790 



Figure 5.— Convergence of the cantilever beam tip displacement. 




(a) Geometry and boundary conditions. 

(b) Finite element model. 

Figure 6.— Clamped rectangular plate analysis— Case II. 

Case II— Analysis of a Clamped Rectangular Plate 

The rectangular plate under a transverse concentrated load, 
shown in figure 6, represents another typical finite element 
test problem. The plate is made of an isotropic material, and 
its parameters are as follows: 


Length, a in 24 

Width, b, in 12 

Thickness, h, in 0.25 

Young’s modulus, E, ksi 30 000 

Poisson’s ratio, v 0.3 

Magnitude of transverse concentrated 
load at center, P, lb 1 000 


The plate is clamped, that is, displacements and rotations are 
restrained, along all four edges. The theoretical solution for 
the transverse displacement at the center (ref. 41) is 

h z — 2.42 X 10~ 3 in. (11) 
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TABLE XI. -COMPARISON OF IFM, MSC/NASTRAN, AND MHOST 
SH75 NORMALIZED CENTER DISPLACMENTS OF 
A CLAMPED RECTANGULAR PLATE 


Number of 
elements, 
(rtXw) 

Normalized displacement 

GIFT/IFM 

PLB4SP 

GIFT/IFMD 

PLB4SP 

MSC/NASTRAN 

QUAD-4 

MHOST 

SH75 

Unite rated 
solution 

Iterated 

solution 

4 (2x2) 

0.825 

0.825 

0.184 

0.006 

0.006 

8 (4x2) 

.987 

.987 

.306 

.010 

.010 

16 (4x4) 

.988 

.988 

.859 

.712 

.726 

32 (8x4) 

1.000 

1. 000 

.945 

.833 

.858 

64 (8x8) 

1 .000 

1.000 

.997 

.953 

.982 


TABLE XII. — MHOST SH75 
NORMALIZED CENTER 
DISPLACMENTS FOR 
SIMPLY SUPPORTED 
RECTANGULAR PLATE 


Number of 
elements, 
(flXttl) 

Normalized 

displacement 

Uniterated 

solution 

Iterated 

solution 

4 (2x2) 

1.279 


8 (4x2) 

.739 

0.818 

16 (4x4) 

.806 

.997 

32 (8x4) 

.799 

.983 

64 (8x8) 

.812 

.913 


The quadrilateral elements of the various software (i.e., 
PLB4SP of GIFT, QUAD-4 of MSC/NASTRAN, QUAD-4 
of ASK A, and SH75 of MHOST) were used to solve the 
problem. In this case, only the center transverse displacements 
are compared for the various methods. The normalized values 
are given in table XI, and the MHOST results obtained for 
the same plate, but with simply supported boundary conditions, 
are given in table XII. 

Case III— Analysis of a Clamped Square Plate by 
Quadrilateral Elements 

A clamped 24-in. square plate, with other parameters 
identical to test Case II, was also analyzed by the quadrilateral 
elements of MSC/NASTRAN, ASKA, GIFT, and MHOST 
as in Case II. The moment resultant at point B (see fig. 6) 
as obtained by the different methods is given in table XIII and 
depicted in figure 7. The transverse center displacements 
computed by the various methods were qualitatively graded 
by the criterion proposed by MacNeal and Harder (ref. 42). 
In their scheme, results are graded as follows on the basis of 
errors in the nodal displacements: 

Grade A less than 2 percent error 

Grade B greater than 2 but less than 10 percent error 


TABLE XIII. -COMPARISON OF IFM, MSC/NASTRAN, AND 
MHOST SH75 NORMALIZED BENDING MOMENTS 


Number of 
elements per 
quarter plate 

Normalized bending moment 

IFM/IFMD 

PLB4SP 

MSC/NASTRAN 

QUAD-4 

MHOST, SH75 

Uniterated 

solution 

Iterated 

solution 

1 

(1 x 1) 

L200 

0 

0 

0 

4 

(2x2) 

.994 

.787 

.620 

1.1137 

9 

(3x3) 

.995 

.875 

.652 

.716 

16 

(4x4) 

.995 

.931 

.732 

1.034 

25 

(5x5) 

— 

— 

.764 

.907 

36 

(6x6) 

— 

— 

.796 

.971 

49 

(7x7) 

— 

— 

.811 

.939 

64 

(8x8) 

— 

— 

.843 

.970 

81 

(9x9) 

— 

— 

.860 

.970 

100 (10x10) 

— 

— 

.860 

.970 

400 (20x20) 

— 

— 

.923 

.986 


Grade C greater than 10 but less than 20 percent error 
Grade D greater than 20 but less than 50 percent error 
Grade F greater than 50 percent error 
For Case III, the grades achieved by the different methods 
are presented in table XIV, and the convergence trend of the 
center transverse displacement with respect to the number of 
elements in the model is depicted in figure 8. 

Case TV— Analysis of a Clamped Square Plate by 
Triangular Elements 

The computations for the clamped square plate of Case III 
were repeated with the triangular plate-bending elements 
PLB3SP of GIFT, TRIA-3 of MSC/NASTRAN, and TRIB-3 
and TUBA-3 of ASKA. The TUBA-3 element of ASKA is 
a higher order element, as described earlier. Results obtained 
from the different methods were again qualitatively graded 
according to the MacNeal and Harder criterion. The grades 
are presented in table XV, and the center transverse 
displacement convergence trend is depicted in figure 9. 
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TABLE XIV. -REPORT CARD FOR QUADILATERAL 
ELEMENTS USED TO SOLVE CLAMPED SQUARE 
PLATE CENTER DISPLACEMENTS 


Number of 
elements 
for full 
plate 
(nXm) 

GIFT/IFM 

GIFT/IFMD 

(PLB4SP) 

MSC/NASTRAN 

(QUAD-4) 

ASKA 

(QUAD-4) 

4 (2x2) 

A 

F 

F 

16 (4x4) 

A 

B 

B 

36 (6x6) 

A 

A 

— 

64 (8x8) 

— 

A 

— 

100 (10x10) 

— 

— 

— 


Method 



Figure 8.— Rate of convergence for rectangular elements used to calculate 
clamped square plate transverse center displacement. 


TABLE XV. -REPORT CARD FOR TRIANGULAR ELEMENTS USED 
TO SOLVE CLAMPED SQUARE PLATE CENTER DISPLACEMENTS 


Number of 
elements 
for full 
plate 

GIFT/IFM 

GIFT/IFMD 

(PLB3SP) 

MSC/NASTRAN 

(TRIA-3) 

ASKA 

(TRIB-3) 

ASKA 

(TUBA-3) 

4 

B 

F 



F 

8 

A 

D 

F 

F 

16 

A 

C 

— 

— 

32 

— 

B 

C 

D 

128 

— 

— 

B 

B 


Method 


Timoshenko 

— O— GIFT/1FM PLB3SP 


— a— ASKATRIB-3 
— A — ASKA TUBA-3 



Number of elements, log (n) 


Figure 9.— Rate of convergence for triangular elements used to calculate 
clamped square plate transverse center displacement. 
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TABLE XVI.— REPORT CARD FOR IFM 
RECTANGULAR ELEMENTS WITH 
DIFFERENT ASPECT RATIOS 


Number of 
elements 
for full 
plate 
(nxm) 

Aspect 

ratio 

Clamped 

boundary 

Simply 

supported 

boundary 

4 (2x2) 

LOO 

A 

A 

4 (2x2) 

1.20 

A 

— 


1.40 

A 

— 


1.60 

A 

— 


1.80 

B 

— 


2.00 

B 

— 

8 (2x4) 

2.00 

A 

— 

16 (4x4) 

LOO 

A 

A 


TABLE XVII, -REPORT CARD FOR HYBRID METHOD 
RECTANGULAR ELEMENTS ON CLAMPED 
SQUARE PLATE 


Number of 
elements 
for full 
plate 

(flXffl) 

GIFT/1FM 

GIFT/IFMD 

PLB4SP 

Mixed method 
MHOST SH75 

Hybrid method 
HMPLS [33] 

4 (2x2) 

A 

F 

F 

16 (4x4) 

A 

C 

C 

64 (8x8) 

A 

B 

A 


Case V— Analysis of Rectangular Plates With Various 
Aspect Ratios 

Rectangular plates identical to Case II, with both simply 
supported and clamped boundary conditions and also with 
different aspect ratios, were examined with the PLB4SP 
element of the GIFT program. Results are given in table XVT. 

Case VI— Analysis of a Clamped Square Plate by 
the Hybrid Method 

This test problem is identical to that of Case IV. The plate 
was examined with the quadrilateral elements of GIFT and 
MHOST and the results compared to the hybrid method 
solution of Chang from element HMPL5 (ref. 33). The 
qualitative grades achieved by the various methods are given 
in table XVII. 


Discussions 

Uniqueness of Elasticity Solution 

In the strict mathematical sense, elasticity solutions are 
unique, that is, for a given force field there is a unique 
displacement state and vice versa. The IFM and IFMD, being 
theoretically equivalent, comply with the uniqueness principle 


(i.e., both IFM and IFMD yield the same solutions; see eq. 5). 
Although initially the IFM and IFMD results are separately 
depicted (tables II to XI), thereafter no distinction is made 
between IFM and IFMD as far as displacement or force 
solutions are concerned. 


Equilibrium Imbalance at the Nodal Points 

To examine the extent to which the different analysis 
methods satisfy the equilibrium conditions at the nodal points 
of a finite element model, the problem of the cantilever beam 
in Case I is considered. The normalized equilibrium imbalance 
at point A (see fig. 4) is defined as 




(F l - F r ) x 100 


( 12 ) 


where F L is the member force at point A from the element 
to the left, F r is the member force at point A from the 
element to the right, and F 0 is the theoretical value for the 
force at point A. 

Tables VI and VII show the error (equilibrium imbalance) 
at point A of the beam model as obtained by GIFT element 
MEMRSP and MHOST elements SH75 and PS 15 1 . The solu- 
tions that were obtained by the integrated force methods do 
not exhibit equilibrium imbalance. The error at point A from 
the MHOST uniterated solution decreases with an increase in 
the number of elements in the discretization. It is about 50 
percent for the model with 2 elements, but about 1 percent 
for a relatively fine model with 24 elements. The results of 
the MHOST iterative scheme, wherein the equilibrium 
imbalance is reduced by a relaxation process, are given in 
tables VI and VII. Note that for this case the MHOST iterated 
scheme has virtually eliminated the error with one iteration. 
The uniterated scheme of MHOST element PS151 does not 
exhibit any error at point A; however, the iterated scheme 
induces minor equilibrium imbalances at that node. 

For the plate flexure problems, the nodal equilibrium 
imbalance is more persistent, especially for the rotational 
degrees of freedom. For Case II, with a 64-element (8x8) 
model, an imbalance of zero is observed in the solution 
obtained by the MSC/NASTRAN QUAD-4 element for the 
transverse translational degree of freedom. However, for the 
rotational degree of freedom the nodal imbalance is of the order 
of magnitude of the reactions developed at the boundary nodes. 
The solution by the IFM (GIFT, element PLB4SP) for this 
case does not exhibit any error for either translational or 
rotational degrees of freedom. 


Attributes of the Stiffness and Pseudo-Stiffness Matrices 

The attributes of the stiffness and IFMD pseudo-stiffness 
matrices, given by equations (6) and (3) respectively, are 
compared for a relatively fine model (i.e., 12 elements) of 
the cantilever beam test problem. Selected global stiffness 
coefficients of the MSC/NASTRAN QUAD-4 element and 
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the pseudo-stiffness coefficients of the GIFT MEMRSP 
element are given in table X. Both the stiffness matrix [K] 
of the MSC/NASTRAN QUAD-4 element with membrane 
response only and the pseudo-stiffness matrix [KJ of the 
GIFT MEMRSP element are symmetrical— of 8x8 dimension. 
Both global stiffness matrices retain similar sign and null 
characteristics. Only the magnitudes of the nonzero coefficients 
differ; that is, with the exception of two elements, the magni- 
tudes of the other 14 elements of the 48th row of the global 
stiffness matrix [K] are higher than those of the pseudo-stiffness 
matrix [KJ. In an overall sense, the stiffness matrix appears to 
be somewhat “stiffer” than the pseudo-stiffness matrix. 

Convergence Trends for Membrane Response 

The normalized tip displacements for the cantilever beam 
of Case I, obtained by GIFT MEMRSP, MSC/NASTRAN 
QUAD-4, and MHOST SH75, are presented in table II. The 
displacements are normalized such that unity represents the 
theoretical solution. For Case I, tip displacement convergence 
is achieved by GIFT MEMRSP for models with four or more 
elements. Both MSC/NASTRAN QUAD-4 and MHOST 
uniterated SH75 converge to approximately 92 percent of the 
theoretical solution. For fewer elements (less than 8 elements 
in the model), the MHOST iterated element-SH75 solution is 
superior to the MSC/NASTRAN and MHOST uniterated 
solutions. However, neither MSC/NASTRAN QUAD-4 nor 
MHOST SH75 uniterated or iterated converge at any closer 
than 92 percent of the closed-form solution, even for a fine, 
48-element model. The displacement convergence trends of 
GIFT MEMRSP and MHOST PS151, given in table III, are 
identical. 

The computed bending stresses at the support of the 
cantilever beam (point A in fig. 4) for different discretizations, 
which were obtained by using the MEMRSP element of GIFT 
and the SH75 and PS151 elements of MHOST, are given in 
tables IV and V. The results are normalized with respect to 
the theoretical bending stress, which is given by 


^theoretical j 04/ 

where M is the bending moment at the support (4800 in. -lb), 
y is the distance from neutral plane (1.0 in.), and I is the 
t(P 

moment of inertia — (=1/6 in. 4 ). 

12 

Both the IFM and IFMD GIFT element MEMRSP yield 
identical results. Furthermore, the stress result converges for 
the first model, which has a single element. The mixed method, 
MHOST, exhibits some error in the computed stress, even for 
fine models; the MHOST uniterated 24-element SH75 model 
has an error of 2.1 percent, which is reduced to 1.2 percent 
by the MHOST iterated solution. The MHOST PS 151 results 
show more rapid convergence, but still require a 16-element 


model to achieve an error of less than 2 percent. The com- 
putational penalty of the MHOST iterated solution (normalized 
to the MHOST uniterated solution) is shown in table VIII for 
element SH75 and in table IX for element PS 151. The addi- 
tional computational time required for the iterated solution is 
one to two times that required for the uniterated solution. 

Convergence Trends for Flexure Response 

The displacements calculated for the clamped rectangular 
plate of Case II are presented in table XI. For this problem, 
GIFT element PLB4SP achieved an accuracy of 98.7 percent 
for a model with 8 elements (4x2). For a coarser model with 
only 4 elements (2x2), the error is about 17.5 percent. The 
solution obtained by MSC/NASTRAN QUAD-4 required 
64 elements to achieve an accuracy of 98 percent. For a 
64-element SH75 model, the MHOST uniterated solution 
achieved an accuracy of about 96 percent; the MHOST iterated 
solution is marginally more accurate. The MHOST results for 
a simply supported rectangular plate, given in table XII, show 
good convergence trends, with minor oscillations for the 
uniterated case. 

For the clamped square plate under transverse concentrated 
load at the center (see table XIV), the ASKA QUAD-4 element 
achieved results no better than a grade of B, even for a model 
with 100 elements. The MSC/NASTRAN QUAD-4 element 
required just 36 elements to produce a “grade A” solution. 
The GIFT PLB4SP-element solution for the most coarse 
model, four elements, achieved a grade of A. 

The convergence trends of the clamped rectangular plate 
bending moment M x at point B (fig. 6), as calculated by GIFT 
PLB4SP, MSC/NASTRAN QUAD-4, and MHOST SH75, 
are presented in table XIII and in figure 7. Note that with the 
GIFT PLB4SP element, convergence for M x at location B 
occurs for the second model, which has four elements per 
quarter plate; the first model, with one element per quarter 
plate, exhibited a bending moment about 20 percent higher 
than Timoshenko’s theoretical solution. The MSC/NASTRAN 
QUAD-4 bending moment shows an error of about 7 percent 
in the fine model with 16 elements per quarter plate. The 
MHOST uniterated solution exhibits about a 7-percent error 
for the model with 400 elements per quarter plate (table XIII), 
but for the same discretization, the MHOST iterated version 
shows a 2-percent error. 

The influence of aspect ratio on the convergence charac- 
teristics of a plate flexure problem was examined by the GIFT 
PLB4SP element. Results, presented in table XVI, show that 
the accuracy decreases as the aspect ratio of the element 
increases from unity (square form); however, PLB4SP retains 
an A grade for the four-element model until the aspect ratio 
reaches 1 .6. For an aspect ratio of 2.0, the whole plate required 
eight elements to secure a grade of A. 

The square plate with clamped boundary was analyzed with 
triangular elements PLB3SP of GIFT, TRIA-3 of MSC/ 
NASTRAN, and TRIB-3 and TUBA-3 of ASKA. Results are 
presented in table XV and figure 9. For element PLB3SP, the 
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result is discernible from the analytical solution for the first 
model, which has four elements in the whole plate; even so, 
the results display engineering accuracy. The next model, with 
eight elements, converges, thereby achieving a grade of A. 
None of the MSC/NASTRAN QUAD-4 nor the ASKA 
TRIB-3 and TUBA-3 results could secure a grade of A, even 
for a finely discretized model with 128 elements, as shown 
in table XVI. 

The displacement convergence characteristics of the IFM 
and the MHOST mixed and hybrid formulations for a clamped 
square plate are given in table XVI. The GIFT PLB4SP 
secured a grade of A for a 4-element model, whereas the 
HMPL5 hybrid element (ref. 33) secured the same grade only 
with a 64-element model. The best grade achieved by the 
MHOST SH75 element for this problem was a B. 

Size of Finite Element Models 

To solve structural mechanics problems, current finite 
element applications employ models with a large number of 
elements and degrees of freedom. Such models are henceforth 
referred to as large models. Although larger models (which 
correspond to smaller finite elements) are presumed to yield 
more accurate solutions, in a strict sense this would be true 
only when element size shrinks to a point, or the displacement 
degrees of freedom are infinite, which is beyond computer 
capability. The question then is, How small should the finite 
elements be in a particular region of a structure in order to 
achieve an acceptable level of accuracy in the prediction of 
stresses and deformations? This question has been researched, 
and techniques such as adaptive mesh refinements have been 
developed. Still, no general answer exists, and mesh refine- 
ment is largely governed by experience and intuition. A related 
issue is the large finite element model (with respect to degrees 
of freedom) whose solution requires thousands of routine 


calculations. Such a model can be handled with intensive 
numerical calculations. This is possible because computation 
has become relatively inexpensive owing to advancements in 
digital computer technology and because accuracy in numerical 
calculations has improved. However, when miniaturization is 
the desirable trend in other disciplines (such as computer 
science, communication engineering, etc.) should large finite 
element models from which solutions are extracted by intensive 
computation be pursued? Perhaps a more appropriate course 
of action would be to search for accurate modeling techniques 
that can generate reliable responses with fewer degrees of 
freedom. The search for such models could be the goal of the 
next generation of finite element technology. 

The issue of model size in finite element calculations is 
explored by taking the examples of the cantilever beam (Case I) 
and plate flexure (Case II) problems. The finite element models 
for the two cases, analyzed by different methods, are depicted 
in figures 10 and 1 1 . For the cantilever beam, only four GIFT 
MEMRSP elements were required to secure a grade of A, 
whereas both MSC/NASTRAN QUAD-4 and MHOST SH75 
elements could secure only of grade of B, even for a fine model 
(see fig. 10). For the plate problem, eight GIFT PLB4SP 
elements were required to achieve a grade of A. To achieve 
the same grade, 64 MSC/NASTRAN QUAD-4 elements were 
needed, whereas 64 MHOST SH75 elements could secure only 
a grade of B. (Note: For both membrane and flexure response, 
the GIFT, MSC/NASTRAN, and MHOST elements are 
equivalent. See NUMERICAL RESULTS.) Overall, the IFM 
required a much smaller model, and the stiffness and mixed 
methods a larger model, to achieve an acceptable level of 
convergence. 

Timoshenko used Ritz’s displacement method to solve a 
plate flexure problem with fixed boundary conditions, 
obtaining accurate solutions with few terms in the series. For 
the square-plate convergence. Case III, the IFMD required 
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Figure 10.— Number of membrane response elements used for various methods. 
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8-Element model 64-Element model 


IFM/IFMD secures grade A MSC/NASTRAN secures grade A 

MHOST secures B 

.Figure 1 1 .—Number of flexural response elements used for various methods. 

four elements with three displacement unknowns. The same 
problem required 36 MSC/NASTRAN QUAD-4 elements, 
which corresponds to 75 displacement unknowns. Likewise, 
64 or more ASKA QUAD-4 or MHOST SH75 elements, or 
64 HMPL5 hybrid elements, all of which correspond to more 
than 100 variables, were required for convergence. For this 
plate problem, only IFMD and Ritz’s convergence 
characteristics are similar; that is, both require a similar 
number of unknowns to achieve convergence. 

All of the finite element analysis methods (the IFM, IFMD, 
stiffness method, hybrid method, and mixed method) are 
approximate formulations. The solutions obtained by these 
methods have to be qualified on the basis of indirect criteria 
such as (1) satisfaction of the equilibrium equations, (2) com- 
pliance of the strain compatibility conditions, and 
(3) elimination of discretization errors by way of the finite 
element model refinements. The IFM attempts to bestow 
balanced emphasis on criteria (1) and (2), and it achieves 
criterion (3) by way of mesh refinement. In other words, all 


three criteria that qualify the solution (equilibrium, 
compatibility, and mesh refinement) are incorporated in the 
IFM; consequently, a converged solution should be accurate 
and reliable. None of the other formulations (stiffness, hybrid, 
and mixed) explicitly impose the strain compatibility condition 
(see figs. 1 and 2); therefore, in a strict sense, there is no 
guarantee that solutions generated by these methods will 
always be correct. 

Concluding Remarks 

Overall, on the basis of the examples analyzed, the following 
conclusions can be drawn: 

1. The integrated force method is superior to the stiffness, 
mixed, and hybrid methods. The latter three methods all 
performed at about the same level. 

2. Most potentials of the stiffness, hybrid, and mixed 
methods have been exploited; these methods probably have 
reached the plateau in their development. The integrated force 
method has now been established and its potential remains to 
be explored. 

3. Since all of the finite element methods are approximate 
in nature, we recommend generating solutions both via the 
integrated force method and the stiffness method and then 
comparing them, rather than qualifying the results by 
successive mesh refinements of any one formulation. 


Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, June 12, 1991 
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Appendix A 

A Plate Flexure Example 


The solution procedure of the integrated force method (IFM) 
is illustrated through the example of a flat cantilever plate in 
flexure (see fig. 12). The plate is made of an isotropic material, 
has a Young’s modulus E of 30 000 ksi, and a Poisson’s ratio 
v of 0.3. The plate is discretized into two rectangular elements, 
each of which has three force and three displacement degrees 
of freedoms per node; the force variables are two moments 
and a shear force, and the displacement variables are two 
rotations and a transverse translation (see fig. 3). 

Solution by the Integrated Force Method 

To analyze the flat cantilever plate by the IFM requires that 
three matrices be generated: the equilibrium matrix [B], the 
flexibility matrix [G], and the compatibility matrix fC] . The 
generation of the matrices is presented in symbolic form to 
avoid algebraic complexity. 


Equilibrium Matrix 

The element equilibrium matrix [B f ,] is the transformation 
that maps nodal loads onto the internal forces at the element 
level. The element equilibrium matrix is a rectangular matrix; 
its rows correspond to the displacement degrees of freedom 
and its columns correspond to independent force variables. 
The consistent equilibrium matrix is generated from the varia- 
tional functional of the IFM. The portion of the functional 
(ref. 5) that yields the matrix [BJ can be written as 



where U pb is the strain energy in flexure; M x , M y , and M xy 



Xi 

- — 



x 2 





Be 

- 

,f 9 . 

X12 

_ _ 




Element equilibrium matrix of 
dimension (12 x 9) 



Boundary compatibility matrix for 
edge AB of dimension (3 x 9) 

(c) 


(a) Displacement degrees of freedom. 

(b) Concatenated displacement. 

(c) Element matrices. 

Figure 12.— Flat cantilever plate in flexture. 
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(A5b) 


d 2 w d 2 w d 2 w 
are the plate-bending moments; and —7, — and — 

, , ® cbr dy oxo y 

are the plate curvatures. J 

The plate domain 0 is defined in a rectangular Cartesian 


coordinate system (x,y). 

The discretized internal energy for the rectangular element 
is expressed as 


*4* = fX] 7 [B,,][F) (A2) 

where C/ dis is discretized internal energy for flexural 
response, [BJ is the plate flexure element equilibrium matrix, 
[X] is the displacement vector of the element, and [F] is the 
force vector of the element. 

The expression to generate a consistent equilibrium matrix 
[BJ is obtained by equating the strain energies given by 
equations (Al) and (A2): 

U ph = t/ d is (A3) 


M v = F 5 4 F b x 4 F-jy + F 8 xy 


M, v = F 9 (A5c) 

The normal moments vary linearly within the element, 
whereas the twisting moment is constant. The constant twisting 
moment M lv will produce interelement discontinuities, which 
of course, if required, can easily be alleviated by a higher order 
polynomial. The assumed moments satisfy the previously 
stated mandatory requirements. 

The element equilibrium matrix is obtained by substituting 
the moments from equation (A5) and the displacements from 
equation (A4) into the energy expression given by equations 
(Al) to (A3) and carrying out the integration. The rectangular 
element equilibrium matrix [BJ is of dimension 12x9; its 
rows correspond to the 12 unknown displacements 

(Xj = W h X 2 - 0vm ^3 = ©w f° r nodes * ~ * t0 4 ) shown in 
figure 12, and its columns correspond to the 9 independent 
force unknowns given by equation (A5). 


The generation of the consistent clement equilibrium matrix 
[BJ requires both the displacement and force distributions in 
the plate domain. For the displacement field, a polynomial 
shape function is chosen in terms of 12 unknowns that satisfy 
the normal plate flexure continuity conditions: 

w(jc f y) = «i 4 ocjX 4 o ?3 y 4 ot 4 x “ 4- cc$xy 4 

+ a 8 .r z y + a 9 xy 2 + a 10 y 3 + «i i* 3 )’ + a l2 xy 3 

(A4) 

The 12 constants (a|,of 2 »--*» a i2) polynomial are 

linked to the 12 nodal displacement degrees of freedom 
(X|,X 2 ,...,X 12 ) of the element by following standard 
techniques. 

Two mandatory requirements of the assumed force field at 
the element level are (1) the force field must satisfy the 
homogeneous equilibrium equation, here, the plate bending 

equation + 2—2 = o) ; and (2) the force 

4 \ dx d >’ dxd >’ / 
components F k (see eq. (A5)) must be independent of one 
another. The latter condition ensures the kinematic stability 
of the element. It is not mandatory that the assumed forces 
satisfy the field compatibility conditions a priori. 

The rectangular element can have 12 nodal forces— 2 
moments and a shear force for each of its 4 nodes. Overall, 
these 12 force components must satisfy the 3 kinematic 
equilibrium conditions; in consequence there are only 9 
independent forces. 

The moment functions of the rectangular element are defined 
in terms of the nine independent force components as 

M x - F, + F 2 x 4 F$y 4 F 4 xy (A5a) 


Flexibility Matrix 

The element flexibility matrix [GJ relates the deformations 
10} to forces [F] as [0] = [GJ[F]. The flexibility matrix is 
symmetrical, of dimension 9x9. Tt is obtained by following 
standard techniques to discretize the complementary strain 
energy U i y which is given as 


U r = ( l/2)[Fj 7 [G,,][Fj = (1/2D) j [m; + M; — 2 vM x M, 


+ ( 1 + v)M~ y 


dx dy 


(A6) 


where D is the flexural rigidity defined as D = (Elrl 12), E 
is Young’s modulus, v is Poisson’s ratio, and h is the plate 
thickness. 

Substituting into equation (A6) the moments M x M y M xy , in 
terms offerees (F„F 2 ,...,F 9 ) as given by equation (A5), and 
integrating yields the 9x9 symmetric flexibility matrix [GJ. 


Compatibility Matrix 

For simplicity, a restrictive procedure to derive the 
compatibility conditions, which is adequate for the plate flexure 
problem, is given here. Generating the compatibility matrix, 
unfortunately, is not as straightforward as generating the 
equilibrium or the flexibility matrices. Refer to references 8, 
10, and 1 1 for the generation of the compatibility conditions 
for finite element analysis. 

The procedure presented here involves direct discretization 
of the continuum plate boundary compatibility conditions by 
using Green’s theorem (ref. 43) and Galerkian’s technique. 
The equation form of the compatibility conditions depends on 
whether such conditions are written for the field or the boundary 
of the elastic domain, since the compatibility principle is unique. 
The field compatibility conditions are incorporated into the 
field integral portion of Green’s theorem, and the correspond- 
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ing boundary compatibility conditions are recovered from the 
boundary integral portion. The Green’s theorem in two 
dimensions can be written as 


ff 

d(j) d\p 

J Jo 

dx dy 


cixdy 



[<j>? + \M dt 


(A7) 


where i and m are the direction cosines to the outward normal 
to the boundary curve. The symbols <f> and \j/ represent 
continuous and differentiable functions of coordinates x and y . 

The plate flexure problem is two-degree indeterminate since 
it has three unknown moments (M x ,M y ,M xy ) but only one field 
equilibrium equation, which is given by 


a 2 M v d 2 M v d 2 M xv 

dx 2 dy 1 dxdy 


(A8a) 


Along the edge where Y = constant, f = 0, and m = 1 , the 
condition given by equation (A1 la) yields two equations: 

(F, + F 3 b) - v(F 3 + F-jb) = 0 (A 12a) 

(F 2 + F 4 b) - v(F 6 + F % b) = 0 (A 12a) 

The condition given by eq. (A1 lb) yields one equation: 

M 9 = 0 (A 12c) 

The three compatibility conditions given by equation (A 12) 
are representative only in the context of discrete finite element 
analysis, because lumped nodal quantities and Galerkian 
integration has not been carried out. The intention here is to 
demonstrate that there are three compatibility conditions per 
edge of the element. The element compatibility condition for 
the edge can be written in matrix form as 


where q is the transverse distributed load. The problem has 
two field compatibility conditions (refs. 5 and 17) given by 


d(M y - vM x ) (1 + v)dM xs _ 
dx dy 

(A8b) 

d(M y - vM y ) (1 +v)dM„_ 
dy dx 

(A8c) 

When Green’s theorem is applied to each of these two field 

compatibility conditions, two boundary 
conditions are recovered: 

compatibility 

(M y - vM x )V - (IF v)M xy m = 0 

(A9a) 

(M x — vM y ) m — ( 1 — v)M X} £ = 0 

(A9b) 

The conditions, specialized for the boundaries of the 
rectangular element, have the following forms: 

I . Along the edges where X = constant, ( = 

1 , and m = 0, 

(M y - vMJ = 0 

(A 10a) 

(M xy ) = 0 

(A 10b) 

2. Along the edges where Y = constant, ( = 

0, and m = 1 , 

(M x - vM x ) = 0 

(Alla) 

(M xy ) = 0 

(A lib) 


The element compatibility conditions in symbolic form are 
obtained by substituting the moment functions (eq. (A5)) into 
the boundary compatibility conditions (eq. (All)). 


[CJ{Fj-[0] (A 13) 

where [CJ is the 3 x9 element boundary compatibility matrix 
for the edge where Y is constant. 

Similar compatibility conditions can be written for the 
element boundary where X is constant, f = 1 , and m — 0. The 
boundary compatibility equation given by equation (A 13) is 
in terms of nine independent forces [Fj; it represents the 
composite compatibility conditions ([C][Gl[F], where 
[CJ = [C][G]), of the IFM for finite element analysis. In the 
computer code GIFT, however, the compatibility matrix [C] 
and the flexibility matrix [G] are generated separately, and 
their product is explicitly determined. The compatibility 
conditions given by equation (A 13), and those obtained by 
generating [C] and [G] separately and taking their product 
[C][G], will have similar characteristics such as bandwidth 
and sparsity, but may be different with respect to some scaling 
factors. 

Integrated Force Method Equations for the Problem 

Each element has 9 independent unknown forces; therefore 

the 2-element discretization has 1 8 force unknowns 1 I , 

mi 

which represent the concatenation of the element forces as 

own /F _ 

j v\ — F]^i,...,F| 2 — F 9e |.\F[o 

= Fid = y (A 14) 

where the subscript iej indicates the /th force of element j. 

The governing equation [S](F] = [PJ for the problem is 
presented next, in equation (A 15). Equation (A 15) contains 
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a total of 18 equations, consisting of 12 equilibrium equations 
(EE) and 6 compatibility conditions (CC), from which the 18 
unknowns can be computed: 



(A 15) 

where 

[B] / equilibrium matrix of dimension 12x9 

for element / 

P\,P 2 ,. .,P \2 12-component mechanical loads 

(5R i ,<5R 2 , 5 R 3 , 

5Rj 6 ,5R 17 ,5R| 8 ) 6-component initial loads 

[C] Mfi edge AB elemental compatibility matrix 

of dimension 3x9 for element / 

The 12 equilibrium equations which link the 18 unknown 

forces to the 12 external loads [P] are assembled 

from the element matrices. The rectangular system of 12 
equilibrium equations has the form 

|[B|1 : [B 2 ]] = [P] (A 16) 

The 12 system equilibrium equations given by equation 
(A 16) occupy the central portion in the IFM governing 
equations given by equation (A 15). 

Because there are 18 force unknowns but only 12 equilib- 
rium equations are available, the plate flexure problem requires 
6 compatibility conditions. These six conditions can be iden- 


tified as the three compatibility constraints along the plate’s 
fixed boundary AB for element 1 and three deformation balance 
conditions for the boundary CD that is common to both 
elements 1 and 2 (fig. 12). 

The three compatibility conditions along boundary AB can 
be written in symbolic form as 

[C,][F,)=0 (A 17) 

The matrix [C]] has a dimension of 3 x9 and is obtained by 
appropriate substitution of direction cosines of equation (A 13) 
for element 1. These three compatibility conditions occupy 
the top position in the IFM governing equation depicted in 
equation (A 15). 

The three compatibility conditions for the common boundary 
CD are given by equation (A1 8). The composite compatibility 
matrix [[C t ] : [C 2 ]] has a dimension of 3 X 1 8 and is obtained 
from element matrices with appropriate assembly for the edge 
CD that is common to elements 1 and 2 (see fig. 12). 

[[C,]:[C 2 ]]j|ijj ={0j (A 18) 

The compatibility condition along interface CD is at the 
bottom location in the IFM governing equation (eq. (A 15)). 
The solution of this governing equation, which contains 
12 EE’s and 6 CC’s, yields the 18 unknown forces. The 12 
displacements can then be obtained from the forces by back 
substitution into equation (2) of the IFM. 

The two-element, finite element solution for the plate flexure 
problem is given in table XVIII along with the strength of 
material beam solutions, which are obtained from a beam 
idealization. Note that the two-element solution yields correct 
moments that are continuous along the interelement boundary 
CD. The maximum transverse displacement obtained for the 
two-element model has only 4.5-percent error compared to 
the theoretical beam solution. 

Solution by the Stiffness Method 

The cantilever plate flexure problem was also solved by the 
stiffness method for the purpose of comparison. The stiffness 
equations are well-known but complicated; therefore, as 
before, the analysis is carried out in symbolic form. To estab- 
lish parallelism between the integrated force and the stiffness 
methods, a slightly different procedure from the normal is 
followed; the purpose will become evident in the process of 
the solution. For the problem, a displacement vector [XJ of 
dimension 24, which represents the concatenation of the 
2 element displacement degrees of freedom, is defined as 

(X f } = (X„ =X um ,.. 1 X, I2 = X I 2h:X ( , 3 

= X I , 2f ...,X, 34 = X l2f2 ) (A 19) 
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TABLE XVIII. -BENDTNG MOMENTS FOR CANTILEVER BEAM 
[See Fig. 12.] 


Nodes 

Deflection, 3 

in. 

M x clement l, a 
in. -lb 

M x element 2, 
in. -lb 

M y element 1 , 
in. -lb 

M y element 2, 
in. -lb 



600.0 

(600.0) 


378.5 


i 




2 


600.0 

(600.0) 


378.5 





3 

.2138 

-300.0 

+ 300.0 

58.60 

-58.60 

4 

.2138 

-300.0 

+ 300.0 

58.60 

-58.60 

5 

.7048 

(.7328) 


0 






6 

.7048 

(.7328) 


0 







Quantities in parentheses arc from the beam solution. 


where the subscript iej represents the /th displacement for the X L ^ — X ]eb — 0 (A2 If) 

jth element (see fig. 12(c)). 

Notice the similarities between the displacement vector [X f .] X i7 = X ]el — X 2e 4 (A2 1 g) 

given by equation (A 19) and the force vector [Fj given by 

equation (A15). These vectors ((Xj and [F]) represent the X c% = X If8 = X 2e $ (A2 1 h) 

concatenation of the element displacement and force degrees 

of freedoms, respectively. By following standard techniques, X& ~ X ]e9 — X le ( 3 (A2 1 i) 

the equilibrium equations given by equation (A 18) can be 

written in terms of nodal displacements [X t .] as X i W = X ]e]0 = X 2e \ (A21j) 


[[Kj] : [K 2 ]]fXJ = [Pj (A20) 


X c \\ = x XcU — x : 


2e2 


(A2 1 k) 


The stiffness matrix fKj] is of dimension 6x12, and its six 
rows represent the contributions to the system equilibrium at 
nodes 3 and 4 (fig, 12). Likewise, the stiffness matrix [K 2 ] 
is of dimension 12x12, and its 12 rows represent the contribu- 
tions to the system equilibrium at nodes 3, 4, 5, and 6. The 
equilibrium equations (eq. (A20)) expressed in terms of 
displacements still represent an indeterminate rectangular 
system with 12 equations in terms of 24 unknown displacement 
variables. Twelve displacement continuity conditions are 
required to augment the equilibrium equations to a solvable 
set of 24 equations in 24 unknowns. The 12 displacement 
continuity conditions for the 2-element plate flexure problem 
are as follows: 


o 

If 

T 

11 

(A21a) 

X* = X U2 = 0 

(A21b) 

X 

II 

>< 

II 

O 

(A21c) 

K-4 — Xif4 — 0 

(A21d) 

+5 = X ie5 = 0 

(A21e) 


X c]2 = X W] 2 ~X M (A211) 

The 12 displacement continuity conditions given by equation 
(A21) can be represented by a single matrix equation: 

[Cry][X c ] = [0] (A22) 

where [C r> ] represents the 12x24 displacement continuity 
matrix. The 12 equilibrium equations (eq. (A20)), written in 
terms of displacements, are coupled to the 12 displacement 
continuity conditions (eq. (A21)) to obtain the 24x24 solvable 
equation system (given by eq. (A23)) of the stiffness method. 
From this system the 24 displacement components {X t .} can 
be calculated: 


[K,]:[K 2 ] 

[Cryl 


IXci = 



(A23) 


The solution of equation (A23), which represents a square 
but nonsymmetrical set of equations, yields the displacements 
from which the forces can be calculated by differentiation, or 
its equivalent, and back calculations. In the popular stiffness 
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method, the continuity conditions (eq. (A2 1)) are trivially 
solved by the linkage of nodal variables and condensation to 
generate the well-known symmetrical stiffness matrix of 
dimension mXm, ( m = 12 for this problem). In principle, 
however, equation (A23) represents the basic unabridged set 
of equations of the stiffness method that, for convenience, is 
manipulated to obtain the condensed symmetrical form. 

From the structure of the IFM equations (eq. (A 15)) and 
the stiffness equations (eq. (A23)), we observe the following: 

(1) In the IFM the equilibrium equations, written in terms 
of forces, are augmented by the compatibility conditions, also 
written in terms of forces, to obtain the IFM governing 
equations [S][F] = [P*], given by equation (A 15). 

(2) In the stiffness method the equilibrium equations, which 
are expressed in terms of displacements, are augmented by 
displacement continuity conditions to obtain the stiffness 
method’s governing matrix equation [K][X] = [P**], given by 
equation (A23). 


(3) For this problem the number of IFM governing equations 
(eq. (A15)) is 18, which is fewer than the 24 equations 
(eq. (A23)) of the displacement method. 

(4) Typically, a sparser system of equations results from 
writing equilibrium equations in terms of forces rather than 
in terms of displacement variables. 

(5) Both the compatibility conditions ([C][G][Fj = [&R]) of 
the IFM and the continuity conditions ([Cry]{X f ] = [0]) of 
the stiffness method yield very sparse systems of equations; 
however, the equations of the continuity conditions are rela- 
tively more sparse than those of the compatibility conditions. 

(6) The equilibrium equations remain indeterminate when 
expressed either in terms of forces or in terms of displacements 
(refer to IFM eq. (A 16) and stiffness eq. (A20)). However, 
the indeterminancy of the equilibrium equations is alleviated 
in the case of the IFM by the compatibility condition 
([C][G](F] = [SR]), or in the displacement method by the dis- 
placement continuity condition [CryHXJ = (0). 
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Appendix B 
Symbols 


[B1 

(m X /?) equilibrium matrix 

U<Us 

discretized internal energy for flexural 

AM 

element equilibrium matrix 


response 

[C] 

(r X n) compatibility matrix 

U P b 

strain energy in flexure 

[CJ 

(3x9) elemental boundary compatibility 

Vy 

shear force at any location along span of beam 


matrix for the edge where y = k\ 

[CJ = [C][G] 

w 

potential of loads 

[C n ] 

displacement continuity matrix 

d 2 w dhv d 2 w 
dx 2 ’ dy 2 ’ dxdy 

plate curvatures 

E 

Young’s modulus 

[Fi 

internal forces; (n x l) internal force vector 

[X! 

nodal displacment unknown 

[FJ 

force component 

[Xr! 

concatenation of 2 elemental nodal 

[G] 

(n x n) concatenated flexibility matrix 


displacements 

[G,.] 

element flexibiity matrix 

Xi,X 2 ...,Xi2 

displacement degrees of freedom 

[G y ] 

(n x m) force coefficient matrix 


constants linked to nodal displacement 

h 

plate thickness 


degrees of freedom 

[j] 

deformation coefficient matrix; first {m x n) 

w 

[G][FJ 


partition of [[S] ] ] T 

[/Jo! 

(n x 1) initial deformation vector 

[K] 

(m x m) symmetrical stiffness matrix 

y.x y y } z>y.xz 

generalized deformation 

[Kv] 

matrix defined by first (w x m) partition of 
[fS][Gl "'[S] 7 ] 

[SR! 

(r x I ) initial deformation vector; 

sr = -\cm 

P,m 

direction cosines to the outward normal to. 

5 V 

displacement at tip of beam 


the boundary curve 

L 

transverse displacement at center of plate 

M x My,M xy 

plate bending moments; generalized stress 
resultants 

[«) 

strain vector 

m 

number of displacement degrees of freedom 


generalized deformations 

[N] 

interpolating polynomials for strain and stress 

K 

curvature 

N x ,N y ,N xy 

generalized stress resultants 

Kz’KyKxy 

generalized deformations 


V 

Poisson’s ratio; 0.3 

n 

number of force degrees of freedom, 



unknown number of equations or forces 

f y 

stress vector 


in IFM 

r xy 

shear stress 

[pi 

{m x 1) external load vector 

<t>,$ 

stress functions in flexure 

[p*i 

equivalent loads 

tz 

transverse rotation 


transverse distributed load 

n 

plate domain in Cartesian coordiantes 

r 

number of compatibility conditions, 

[0] 

null matrix 

[SI 

r — ft — fii 

(nxn) IFM governing matrix 

Superscript: 


u r 

complementary strain energy 

T 

transpose of matrix or vector 
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